#include <arm_neon.h>
#include <math.h>

#define TAYLOR_PRECISION 1e-6

// 计算log函数
double neon_log(double x) {
    if (x <= 0.0) {
        return NAN; // log(0)和log(负数)是未定义的
    }
    float64x1_t one = vdup_n_f64(1.0);
    float64x1_t y = vsub_f64(vdup_n_f64(x), one);
    float64x1_t term = y;
    float64x1_t result = y;
    float64x1_t n = vdup_n_f64(1.0);

    while (vget_lane_f64(vabs_f64(term), 0) > TAYLOR_PRECISION) {
        n = vadd_f64(n, one);
        term = vmul_f64(term, vneg_f64(y));
        result = vadd_f64(result, vdiv_f64(term, n));
    }

    return vget_lane_f64(result, 0);
}

// 计算exp函数
double neon_exp(double x) {
    float64x1_t one = vdup_n_f64(1.0);
    float64x1_t term = one;
    float64x1_t result = one;
    float64x1_t n = vdup_n_f64(1.0);

    while (vget_lane_f64(vabs_f64(term), 0) > TAYLOR_PRECISION) {
        n = vadd_f64(n, one);
        term = vmul_f64(term, vdup_n_f64(x));
        term = vdiv_f64(term, n);
        result = vadd_f64(result, term);
    }

    return vget_lane_f64(result, 0);
}

// 计算pow函数
double neon_pow(double base, double exponent) {
    if (base == 0.0 && exponent <= 0.0) {
        return NAN; // 0的负数次幂和0的0次幂是未定义的
    }
    double log_base = neon_log(base);
    double exp_result = neon_exp(log_base * exponent);
    return exp_result;
}

int main() {
    double bases[] = {2.0, -2.0, 0.0, 2.5, -2.5};
    double exponents[] = {3.0, -3.0, 0.0, 3.7, -3.7};

    for (int i = 0; i < 5; i++) {
        for (int j = 0; j < 5; j++) {
            double base = bases[i];
            double exponent = exponents[j];
            double result = neon_pow(base, exponent);
            printf("pow(%f, %f) = %f\n", base, exponent, result);
        }
    }

    return 0;
}


#include <arm_neon.h>
#include <math.h>

#define TAYLOR_PRECISION 1e-6

// 计算log函数
float neon_log(float x) {
    if (x <= 0.0f) {
        return NAN; // log(0)和log(负数)是未定义的
    }
    float32x2_t one = vdup_n_f32(1.0f);
    float32x2_t y = vsub_f32(vdup_n_f32(x), one);
    float32x2_t term = y;
    float32x2_t result = y;
    float32x2_t n = vdup_n_f32(1.0f);

    while (vmaxv_f32(vabs_f32(term)) > TAYLOR_PRECISION) {
        n = vadd_f32(n, one);
        term = vmul_f32(term, vneg_f32(y));
        result = vadd_f32(result, vdiv_f32(term, n));
    }

    return vget_lane_f32(result, 0);
}

// 计算exp函数
float neon_exp(float x) {
    float32x2_t one = vdup_n_f32(1.0f);
    float32x2_t term = one;
    float32x2_t result = one;
    float32x2_t n = vdup_n_f32(1.0f);

    while (vmaxv_f32(vabs_f32(term)) > TAYLOR_PRECISION) {
        n = vadd_f32(n, one);
        term = vmul_f32(term, vdup_n_f32(x));
        term = vdiv_f32(term, n);
        result = vadd_f32(result, term);
    }

    return vget_lane_f32(result, 0);
}

// 计算pow函数
float neon_pow(float base, float exponent) {
    if (base == 0.0f && exponent <= 0.0f) {
        return NAN; // 0的负数次幂和0的0次幂是未定义的
    }
    float log_base = neon_log(base);
    float exp_result = neon_exp(log_base * exponent);
    return exp_result;
}

int main() {
    float bases[] = {2.0f, -2.0f, 0.0f, 2.5f, -2.5f};
    float exponents[] = {3.0f, -3.0f, 0.0f, 3.7f, -3.7f};

    for (int i = 0; i < 5; i++) {
        for (int j = 0; j < 5; j++) {
            float base = bases[i];
            float exponent = exponents[j];
            float result = neon_pow(base, exponent);
            printf("pow(%f, %f) = %f\n", base, exponent, result);
        }
    }

    return 0;
}
